#ifndef MYCUDA_H
#define MYCUDA_H
#include <stdio.h>
#include <cstdlib>
#include <cuda_runtime.h>
#include <cuda.h>
#include "timer.h"
#include <device_launch_parameters.h>
#include <helper_functions.h>
#define SIZE_X 100
#define SIZE_Y 100
#define SIZE_Z 100
#define _ANGLE 1.57

__device__ static float3 _getPos(int pos)
{
	float3 p;
	int slice=SIZE_X*SIZE_Y;
	p.z=pos/slice;//slice
	pos%=slice;
	p.y=pos/SIZE_X;
	p.x=pos%SIZE_X;
	return p;
}
__device__ static float _getAngle(float3 v21,float3 v23)//get angle between v21 and v23
{
	float dot=v21.x*v23.x+v21.y*v23.y+v21.z*v23.z;
	float d21=sqrt(v21.x*v21.x+v21.y*v21.y+v21.z*v21.z);
	float d23=sqrt(v23.x*v23.x+v23.y*v23.y+v23.z*v23.z);
	return dot/(d21*d23);
}
__device__ static float _compDist(float3 pc,float3 pa,float3 pb)
{
	float SAB=sqrt((pa.x-pb.x)*(pa.x-pb.x)+(pa.y-pb.y)*(pa.y-pb.y)+(pa.z-pb.z)*(pa.z-pb.z));
	float SAC=sqrt((pa.x-pc.x)*(pa.x-pc.x)+(pa.y-pc.y)*(pa.y-pc.y)+(pa.z-pc.z)*(pa.z-pc.z));
	float SBC=sqrt((pc.x-pb.x)*(pc.x-pb.x)+(pc.y-pb.y)*(pc.y-pb.y)+(pc.z-pb.z)*(pc.z-pb.z)  );
	float P=(SAB+SAC+SBC)/2;
	float SS=sqrt(P*(P-SBC)*(P-SAC)*(P-SAB));

	float3 v21,v23;
	v21.x=pc.x-pa.x;	v21.y=pc.y-pa.y;	v21.z=pc.z-pa.z;
	v23.x=pb.x-pa.x;	v23.y=pb.y-pa.y;	v23.z=pb.z-pa.z;

	float angel=_getAngle(v21,v23);
	if(angel<cos(_ANGLE))//angel cab >90,then dist=ca
		return sqrt(v21.x*v21.x+v21.y*v21.y+v21.z*v21.z);
	else 
	{
		v21.x=pc.x-pb.x;	v21.y=pc.y-pb.y;	v21.z=pc.z-pb.z;
		v23.x=pa.x-pb.x;	v23.y=pa.y-pb.y;	v23.z=pa.z-pb.z;
		angel=_getAngle(v21,v23);
		if(angel<cos(_ANGLE))	//angel cba >90,then dist=cb
			return sqrt(v21.x*v21.x+v21.y*v21.y+v21.z*v21.z);
	}
	return (SS+SS)/SAB;
}
__global__ static void sumOfDataTerm(float4 *array1,float4* array2,float* result)
{
    const int tid = threadIdx.x;
    //const int bid = blockIdx.x;

	float4 m1=array1[tid];
	float4 m2=array2[tid];

	float3 mypos=_getPos((int)m1.w);	
	float theta=m2.w;
	float3 pa,pb;
	pa.x=m1.x;	pa.y=m1.y;	pa.z=m1.z;
	pb.x=m2.x;	pb.y=m2.y;	pb.z=m2.z;

	float dis=_compDist(mypos,pa,pb);
	dis*=dis;

	float energy=log(sqrt(2*3.14f)*theta) +dis/(theta*theta);
	if(energy<-500)	energy=-500;
	if(energy>500)	energy=500;

	result[tid]=energy;
}
#endif
